library(specr)
library(fixest)
library(readr)
library(broom.mixed)
library(texreg)
library(dotwhisker)
library(dplyr)
library(fixest)

# Two-Way linear Fixed Effect Models that only use data
# where there is variation
# across time in the outcomes - 'mimicking' logit panel
# regression outcomes

source("funs_and_constants.R")

# add in data
epr_yearly_group_clusters_h1 <- read_csv("data/epro_data_h1.csv")
epro_year_group_clusters <- read_csv("data/epro_data_h2.csv")
epro_year_group_clusters_geo <- read_csv("data/epro_geo_data_h3.csv")

# excluding obs without variation
# strategy: if unique values between country groups is 1, nothing changes
# if there is more than one unique values, this indicates change within
# that particular country

epr_yearly_group_clusters_h1 <- epr_yearly_group_clusters_h1 %>%
    group_by(countries_gwid) %>%
    # also filter nas because missing rows can cause rows with variation to not be included
    filter(if_any(c(disagreement, share_disagreement, n_orgs, resource_and_agriculture, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, n_tek_groups), is.na)==FALSE) %>%
    mutate(n_values_agr = length(unique(disagreement)),
           variation_disagreement = if_else(n_values_agr > 1, TRUE, FALSE),
           n_values_prop = length(unique(share_disagreement)),
           variation_share_disagreement = if_else(n_values_prop > 1, TRUE, FALSE))

epro_year_group_clusters <- epro_year_group_clusters %>%
  group_by(countries_gwid) %>%
  # also filter nas because missing rows can cause rows with variation to not be included
  filter(if_any(c(disagreement, share_disagreement, n_orgs, religious_segments_bin, n_ed_religions, hhi_rel, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic, n_tek_groups), is.na)==FALSE) %>%
  mutate(n_values_agr = length(unique(disagreement)),
         variation_disagreement = if_else(n_values_agr > 1, TRUE, FALSE),
         n_values_prop = length(unique(share_disagreement)),
         variation_share_disagreement = if_else(n_values_prop > 1, TRUE, FALSE))


epro_year_group_clusters_geo <- epro_year_group_clusters_geo %>%
  group_by(countries_gwid) %>%
  # also filter nas because missing rows can cause rows with variation to not be included
  filter(if_any(c(disagreement, share_disagreement, n_orgs, n_non_intersection_group_polygons,
                  multiple_polygons_bin, nightlight_total_pc_log, n_tek_groups,
                  spatial_hhi, regaut, groupsize, incidence_flag, status_pwrrank, any_multiethnic), is.na)==FALSE) %>%
  mutate(n_values_agr = length(unique(disagreement)),
         variation_disagreement = if_else(n_values_agr > 1, TRUE, FALSE),
         n_values_prop = length(unique(share_disagreement)),
         variation_share_disagreement = if_else(n_values_prop > 1, TRUE, FALSE))

# to compare - n obs should be exactly the same
# full_model_h1_log <- feglm(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
#                        | year + countries_gwid, cluster = ~ countries_gwid,
#                        data = epr_yearly_group_clusters_h1, family = "binomial")


full_model_h1 <- feols(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                       data = epr_yearly_group_clusters_h1 %>% filter(variation_disagreement == TRUE))

full_model_h1 <- feols(disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year + countries_gwid, cluster = ~ countries_gwid,
                       data = epr_yearly_group_clusters_h1)

full_model_h1_max <- feols(share_disagreement ~ resource_and_agriculture + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic +  n_tek_groups |
                             year + countries_gwid, cluster = ~countries_gwid,
                           data = epr_yearly_group_clusters_h1 %>% filter(variation_share_disagreement == TRUE))

full_model_h1_n_resources <- feols(disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                                   +  n_tek_groups | year+ countries_gwid, cluster = ~countries_gwid,
                                   data = epr_yearly_group_clusters_h1 %>% filter(variation_disagreement == TRUE))

full_model_h1_n_resources_max <- feols(share_disagreement ~ none_one_or_both_nats_agri_char + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic +  n_tek_groups | 
                                         year + countries_gwid,
                                       data = epr_yearly_group_clusters_h1 %>% filter(variation_share_disagreement == TRUE),
                                       cluster = ~countries_gwid)

model_list_h1 <- list(full_model_h1, full_model_h1_max, full_model_h1_n_resources, full_model_h1_n_resources_max)

texreg::screenreg(model_list_h1, stars = stars, custom.coef.map = coef_map)

# old table
# texreg(model_list_h1,
#        stars = stars,
#        custom.coef.map = coef_map,
#        table = FALSE,
#        custom.model.names = c("Disagreement (1/0)", "Prop. Disagreement", "Disagreement (1/0)", "Prop. Disagreement"),
#        include.proj.stats = FALSE,
#        file = "tables/h1_lpm_panel_no_all_zero.tex")

texreg(list(full_model_h1, full_model_h1_n_resources, full_model_h1_max, full_model_h1_n_resources_max),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:2, "Prop. Disagreement" = 3:4),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h1_lpm_panel_no_all_zero.tex")


####
## h2
#####

full_model_h2 <- feols(disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic +
                        n_tek_groups | year +countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(variation_disagreement == TRUE))

# n rels as IV
full_model_h2_n_rels <- feols(disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                              +  n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(variation_disagreement == TRUE))


full_model_h2_herf <- feols(disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                            +  n_tek_groups | year + countries_gwid, cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(variation_disagreement == TRUE))

full_model_h2_max <- feols(share_disagreement ~ religious_segments_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic +  n_tek_groups | year + countries_gwid,
                           cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(variation_share_disagreement == TRUE))

full_model_h2_max_n_rels <- feols(share_disagreement ~ n_ed_religions + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic +  n_tek_groups | year + countries_gwid,
                                  cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(variation_share_disagreement == TRUE))

full_model_h2_max_herf <- feols(share_disagreement ~ hhi_rel + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic +  n_tek_groups | year + countries_gwid,
                                cluster = ~countries_gwid, data = epro_year_group_clusters %>% filter(variation_share_disagreement == TRUE))


screenreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf),
          stars = stars)

texreg(list(full_model_h2, full_model_h2_n_rels, full_model_h2_herf, full_model_h2_max, full_model_h2_max_n_rels, full_model_h2_max_herf),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h2_lpm_panel_no_all_zero.tex")


######
### h3
#######

full_model_h3 <- feols(disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                       +  n_tek_groups | year + countries_gwid,
                       cluster = ~countries_gwid,
                       data = epro_year_group_clusters_geo %>% filter(variation_disagreement == TRUE),
)

full_model_h3_max <- feols(share_disagreement ~ n_non_intersection_group_polygons + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                           +  n_tek_groups | year + countries_gwid,
                           cluster = ~countries_gwid,
                           data = epro_year_group_clusters_geo %>% filter(variation_share_disagreement == TRUE),
)

full_model_h3_bin <- feols(disagreement ~ multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                           +  n_tek_groups | year + countries_gwid,
                           cluster = ~countries_gwid,
                           data = epro_year_group_clusters_geo %>% filter(variation_disagreement == TRUE),
)

full_model_h3_max_bin <- feols(share_disagreement ~  multiple_polygons_bin + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                               +  n_tek_groups | year + countries_gwid,
                               cluster = ~countries_gwid,
                               data = epro_year_group_clusters_geo %>% filter(variation_share_disagreement == TRUE),
)

full_model_h3_hhi <- feols(disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                           +  n_tek_groups | year + countries_gwid,
                           cluster = ~countries_gwid,
                           data = epro_year_group_clusters_geo %>% filter(variation_disagreement == TRUE),
)

full_model_h3_max_hhi <- feols(share_disagreement ~ spatial_hhi + n_orgs + regaut + groupsize + incidence_flag + status_pwrrank + any_multiethnic
                               +  n_tek_groups | year +countries_gwid,
                               cluster = ~countries_gwid,
                               data = epro_year_group_clusters_geo %>% filter(variation_share_disagreement == TRUE),
)


screenreg(list(full_model_h3, full_model_h3_max, full_model_h3_hhi, full_model_h3_max_hhi, full_model_h3_bin, full_model_h3_max_bin), stars = stars)

texreg(list(full_model_h3, full_model_h3_bin, full_model_h3_hhi, full_model_h3_max, full_model_h3_max_bin, full_model_h3_max_hhi),
       stars = stars,
       custom.coef.map = coef_map,
       custom.header = list("Disagreement (1/0)" = 1:3, "Prop. Disagreement" = 4:6),
       table = FALSE,
       include.proj.stats = FALSE,
       file = "tables/h3_lpm_panel_no_all_zero.tex")






